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Abstract — We study synchronization of Kuramoto os- 
cillators in strongly modular networks in which the struc- 
ture of the network inside each community is averaged. 
We find that the dynamics of the interacting communities 
can be described as an ensemble of coupled planar oscil- 
lators. In the limit of a large number of communities, we 
find a low dimensional description of the level of synchro- 
nization between the communities. In this limit, we de- 
scribe bifurcations between incoherence, local synchrony, 
and global synchrony. We compare the predictions of this 
simplified model with simulations of heterogeneous net- 
works in which the internal structure of each community 
is preserved and find excellent agreement. Finally, we in- 
vestigate synchronization in networks where several layers 
of communities within communities may be present. 

1. Introduction 

Large networks of coupled oscillators are found in 
many applications of science and engineering, including 
synchronized flashing of fireflies JTJ, cardiac pacemaker 
cells |2, oscillations of pedestrian bridges Q, and circa- 
dian rhythms in mammals (4). The Kuramoto model Q 
has become a paradigm for modeling and studying emer- 
gence of collective behavior in the form of synchro- 
nization. In the Kuramoto model each oscillator is de- 
scribed by a phase angle 0„ that evolves as 6 n = to,, + 
K 2m = i A, ml sm(6 m —6„), where u>„ is the intrinsic frequency 
of oscillator n, K is the global coupling strength, A nm en- 
codes the network topology, and n,m = 1 , . . . , N. 

The path to synchrony in typical non-modular networks 
is characterized by incoherence for small K, followed by 
the emergence of a single synchronized cluster when K sur- 
passes a critical value K c J5J. However, when the network 
structure is modular, synchrony can occur hierarchically: 
first locally within each community, and then globally as 
communities synchronize with one another. Several studies 
on synchronization in modular networks exist |]7] |8] I91H01. 
but few analytic results for large networks of heterogeneous 
oscillators with many communities exist. 

In this paper we study a system of equations previously 
explored in Refs. JU [10] QT) using the dimensionality re- 
duction techniques of Ott and Antonsen ifTTI . We find an- 
alytical expressions for local and global order parameters 
describing synchronization within communities and on the 
whole network, respectively, and completely characterize 



the phase space of the system. 

This paper is organized as follows. In Sec. 2 we sum- 
marize the results from Ref. [1101 where network structure 
within communities is averaged. In Sec. 3 we compare 
the results of the averaged system to numerical simulations 
where the structure of each community is preserved. In 
Sec. 4 we show that results from Sec. 2 generalize to sys- 
tems with multiple levels of community structure by study- 
ing the case of three levels. In Sec. 5 we conclude with a 
brief discussion. 

2. Hierarchical Synchrony in Two-level Hierarchical 
Networks 

We consider a system of C communities, labeled cr - 
1, . . . , C, each containing oscillators. Oscillator « in 
community cr has phase ff^ and evolves according to ll8l [TTI 

C i N i 

% = < + V Va . — V sin(< - <£), ( 1 ) 

cn = l m=l 

where co% is its intrinsic frequency, K' Tcr is the coupling 
between oscillators in communities <x and cr', and Jja- is the 
fraction of oscillators in community cr. 

For analytic tractability, we make the following sim- 
plifications. First, we assume all communities are of the 
same size, i.e. N^- = N and rju- = C l . In addition, 
we assume that (i) the coupling strength between oscilla- 
tors within the same community is much larger than the 
coupling strength between oscillators in different commu- 
nities and (ii) the intrinsic frequency for an oscillator is 
drawn from a distribution specific to its own community. 
To ensure condition (i) for a large number of communi- 
ties, we let K trcr — Ck if cr = o J , and K otherwise, where 
k and K are of the same order. Finally, to ensure condi- 
tion (ii) we assume that the frequency to^ is drawn from 
a distribution g a {co). We will assume that this distribu- 
tion is Lorentzian with uniform spread 6 and community- 
specific mean Q. a : go-(w) = 7t~ l 6/[6 2 + (co — Q. (r ) 2 '\. Further- 
more, the means Q. a are drawn from their own distribution 
G{Q) which is Lorentzian with spread A and mean zero: 
G(Q) = 7r*A/(A 2 + Q 2 ). Parameters 6, A = 1 are used for 
all figures presented. A discussion of how some of these 
assumptions can be relaxed can be found in Ref. IfTOl . 

Finally, to characterize the degree of local and global 
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Figure 1 : Phase space for Eq. ([TJ with 5, A = 1 . Regions 
A, B, C, and D (described in the text) are denoted in red, 
yellow, green, and blue, respectively, with bifurcations (i)- 
(iv) indicated by solid and dashed curves. 

synchrony, we define the complex order parameters 

N C 

^ = r ^i|y», z=^ = i^ Z(rj ( 2) 
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In regions A and B, where R — 0, the communities de- 
couple and the steady-state degree of local synchrony for 
each community is given by 
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k-K/C 



if k-K/C < 25, 
otherwise, 



(5) 



and the onset of local synchrony is given by bifurcation (i) 
k - K/C = 25. 

Now we analyze regions C and D. Given a steady-state 
R value (to be discussed later) it can be shown [10| that 
in region C all communities lock and their degree of local 
synchrony is given implicitly by 
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(6) 



In region D only a fraction of the the communities lock, 
which are precisely those communities with sufficiently 
small mean frequency (in magnitude) given by 



IQJ <D. = KR 



such that r a measures the degree of local synchrony in 
community <x and R measures the degree of global syn- 
chrony over the entire network. To measure the average 
degree of local synchrony we introduce 7 — ^ £o- r LT . 

The phase space of this system as a function of K, k was 
found in IflOl in the limits N, C — > oo and is shown in Fig- 
ureQ] For small K, k (A: red) both the average local degree 
of synchrony 7 and global degree of synchrony R are zero. 
If k is increased and K is kept sufficiently small (B: yellow) 
local synchrony is non-zero but no global synchrony exists. 
Finally, for sufficiently large K (C: green and D: blue) both 
local and global synchrony are non-zero. As we will dis- 
cuss, these regions differ in how communities synchronize 
with one another. Bifurcation curves between these differ- 
ent states are labelled (i)-(iv) and will be discussed below. 
We now summarize some results from Ref. IfTOl on local 
and global synchrony. 

2.1. Local synchrony 

To classify local synchrony we consider the continuum 
limit N — > oo and introduce the disribution fa-(&, <^>, t) that 
describes the density of oscillators with phase and fre- 
quency 9 and to at time t. After a dimensionality reduction 
(described in detail in Refs. IflOl [TP '), we find that r a and 
iff a- evolve according to the 2C-dimensional planar oscilla- 
tor system 



1 



-1/2 



(7) 



and have a degree of local synchrony given by Eq. (|6}. The 
drifting communities turn out to have solutions where 
oscillates with approximate mean VI ~ 2S/(k - K/C). Bi- 
furcation (ii) is given by k-K/C = 25, which can be found 
by sending Q. — > oo in Eq. (0. 

2.2. Global synchrony 

To classify global synchrony we consider the continuum 
limit C — > oo and introduce the distribution F(ifr, Q., r, t) 
that describes the density of communities with phase, mean 
frequency, and degree of local synchrony iff, Q., and r at 
time t. After a dimensionality reduction (described in detail 
in Ref. IflOl ). we find that R evolves according to 
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where r solves Eq. (O for O = -zA, while *F = 0. 
steady-state degree of global synchrony is given by 
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The 
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which can be solved consistently with f. Sending r — > 
y/4-A/K - 1 + yields an onset of global synchrony given by 
k = ~ 7 corresponding to bifurcations (iii) and (iv). 

We now consider the implications of our analysis on the 
hierarchy of local and global synchrony. Starting in region 
A, the onset of local and global synchrony will occur sep- 
arately if bifurcation (i) is reached first, otherwise bifur- 
cation (iii) is reached first and they occur simultaneously. 




Figure 2: Degree of global and local synchrony (blue cir- 
cles and red triangle) along paths (a) k = 3K/2 and (b) 
k — K/2 from simulation with N,C - 1000. Theoretical 
predictions are plotted in dashed black. 

Given a path k = mK where K is increased from zero, if 
m > m c — ; = = then the relative ratio of local 

A-tf+ VA 2 +(5 2 +6(5A 

to global coupling is strong enough to produce hierarchical 
synchrony, where 7 > before R > 0. On the other hand, 
if m < m c global effects dominate at onset and 7 > 0, R > 
occur simultaneously. We plot the steady-state values of 
R and 7 from simulations (blue circles and red triangles) 
and theory (dashed black curves) resulting from increasing 
K while k = mK for N,C = 1000 in Fig. for m = 3/2 
and 1/2 [subfigures (a) and (b), respectively]. We observe 
a hierarchical separation of the onset of local and global 
synchrony for m — 3/2 but not for m = 1/2. 

3. Numerical Experiments 

In this Section, we test the usefulness of the theoreti- 
cal results in Section 2 in predicting behavior of non-trivial 
networks. We consider a system given by 

K c N 

% = < + n Z Z A ™ s[n(e <" - 6n) ' (10) 

o"'=l m=l 

where A™ encodes the network structure between oscilla- 
tors n and m in communities cr and cr', respectively, and K 
is the global coupling strength. We consider the case where 
the network in community cr is an Erdos-Renyi random net- 
work. In order to attain modularity, we assume that when 
cr = cr', A™ is 1 with probability p\ ~ 1 and otherwise, 
and that when cr ± cr', A^ is 1 with probability p 2 <k 1 
and otherwise. Using the mean degrees {d)\o = Np\,2 we 
can estimate effective k,K values as k = {d)\K/N = Kp\ 
and K = C(d) 2 K/N = CKp 2 . 

In Fig. [3] we plot the degree of global and local syn- 
chrony from simulating Eq. (fTUl i on Erdos-Renyi networks 
(blue circles and red triangles) compared to theoretical pre- 
dictions from the community-averaged case (dashed black) 
using the corresponding values of k and K estimated in the 
previous paragraph. Parameters are N = 200, C = 50, 
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Figure 3: Degree of global and local synchrony (blue cir- 
cles and red triangle) for Erdos-Renyi network simulations 
with N = 200, C = 50, pi = 0.95, and p 2 = 0.0126 (a) 
and 0.038 (b), respectively. Theoretical predictions from 
the averaged system are plotted in dashed black. 

pi = 0.95, and p 2 = 0.0126 and 0.038 [subplots (a) and 
(b), respectively]. The agreement between the theoretical 
prediction from the community-averaged system and ac- 
tual simulation of the Erdos-Renyi networks is excellent. 

4. Hierarchical Synchrony in Multiple-Levels 

We now consider a system with three levels of commu- 
nity structure: a network of oscillators with C communities 
where each community contains S sub-communities, and 
each sub-community contains N oscillators that evolve ac- 
cording to 

I c i s K cra ~' N 

@n,y,(T ~ ^n,y,a ^ ^ ^ ~~ ^ ^ ~jy ^ ^ sin($ m y' iCr ' @n,y,cr)> 
cr'=l y'=l m=\ 

where 8 n ^ tCr and a>„ )7jtr denote the phase and frequency of 
oscillator n in sub-community y of community cr and KZy 
denotes the coupling strength between oscillators in sub- 
community y of community cr and sub-community y' in 
community cr' . To ensure multi-level community structure 
we assume that K™, = SCK\ for y — y' ,cr — cr' , CK 2 
for y + y 1 ' ,cr — cr' , and K3 otherwise. In analogy with 
the two-level system, the frequencies Oi n ,y,<r are drawn from 
the Lorentzian f y (cS) = n^d/ld 2 + (oj - Q^) 2 ], is 
drawn from the Lorentzian Go-(Q) = Tr'A/l^-KQ-Qj) 2 ], 
and Q% is drawn from the Lorentzian G(£X) = n~ 1 A/(A 2 + 
Q 2 ). To measure the degrees of local, mid-level, and global 
synchrony, we use the order parameters 

n= 1 y=l o~= 1 

whose magnitudes rj" 7 , r%, and measure the degrees of 
synchrony and have phases ifr Y l °', if/^, and 1^3. We measure 
the average local and mid-level degrees of synchrony by 
h = jc Z y ,o- r T and n = r 2 - 





Figure 4: Degree of global, mid-level, and local synchrony (blue circles, red triangle, and green crosses) for three com- 
munity levels with N = 500, S = 200, C = 120. Coupling is described by paths K2 = aK^, K\ - bKi for a,b — |, | (a), 
|,i(b), I, |(c), and 1.1(d). 



We investigate the hierarchy of local, mid-level, and 
global synchrony by moving through paths K2 = aK^, 
K\ = bK2 as we increase ^3 from zero. In Fig. |4]we plot 
the resulting r$, ?2, and T\ (blue circles, red triangles, and 
green crosses, respectively) for a, b = |, | (a), |, | (b), 
i, I (c), and |, i (d). Onsets of local, mid-level, and global 
synchrony, denoted k\, hi, and £3, respectively, are indi- 
cated by arrows. For a,b = §, | the onsets of each level 
of synchrony are separated, while for a, b = h, j all occur 
simultaneously. For a,b — |, i mid-level and global syn- 
chrony occur simultaneously and are separated from local 
onset while for a,b = \, | local and mid-level synchrony 
occur simultaneously and are separated from global onset. 

5. Discussion 

We have presented analytical results describing local and 
global synchrony in modular networks where the network 
structure within each community is averaged. Furthermore, 
we have shown via numerical simulations that these analyt- 
ical results predict very well the dynamics of non-trivial 
networks with community-wise Erdos-Renyi topologies. 
Importantly, analytical results indicate whether the path 
to synchrony occurs hierarchically or not. The effect of 
stronger heterogeneity on hierarchical synchrony in modu- 
lar networks remains an open area of research. Finally, we 
have investigated hierarchical synchrony in networks with 
several layers of community structure and have found that 
depending on the relative ratios of the coupling strengths, 
synchrony of different layers can occur either hierarchically 
or simultaneously as in the two-layer case. 
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